Hyperbolic metamaterials: nonlocal response regularizes broadband super-singularity 
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We study metamaterials known as hyperbolic media that in the usual local-response approxima- 
tion exhibit hyperbolic dispersion and an associated broadband singularity in the density of states. 
Instead, from the more microscopic hydrodynamic Drude theory we derive qualitatively different 
optical properties of these metamaterials, due to the free-electron nonlocal optical response of their 
metal constituents. We demonstrate that nonlocal response gives rise to a large-wavevector cutoff 
in the dispersion that is inversely proportional to the Fermi velocity of the electron gas, but also for 
small wavevectors we find differences for the hyperbolic dispersion. Moreover, the size of the unit 
cell influences effective parameters of the metamaterial even in the deep sub-wavelength regime. 
Finally, instead of the broadband super-singularity in the local density of states, we predict a large 
but finite maximal enhancement proportional to the inverse cube of the Fermi velocity. 

PACS numbers: 42.70.Qs, 78.20.Ci, 71.45.Gm, 71.45.Lr 



I. INTRODUCTION 

Metamaterials, consisting of subwavelength artificial 
unit cells, show a great potential in optical applications, 
such as perfect lenses^ and invisibility cloaks. 2,3 Hyper- 
bolic metamaterials (HMM) are of special interest be- 
cause of their unusual hyperbolic dispersion curves^— 
that support radiative modes with unbounded wavenum- 
bers. Because of the diverging radiative local density of 
states (LDOS), a point emitter in such a medium would 
exhibit instantaneous radiative decay£~— This broad- 
band 'super-singularity'^ is indeed broadband, since 
the hyperbolic dispersion does not rely on specific res- 
onances. The prediction by Jacob et al& that hyper- 
bolic media thereby form a new route to enhanced light- 
matter coupling recently found experimental support 
These measured lifetimes were nonzero, and state-of-the- 
art theories explain this from three parameters: the non- 
vanishing damping 7, the size a of the unit cell, and 
the size D of the emitter In particular, as a func- 
tion of these parameters the radiative LDOS scales as 
7 - 3 / 2 (Ref. [1), a" 3 (Ref. 7) and D~ 3 (Ref. i), respec- 
tively. Usually unit cells are larger than emitter sizes, 
which makes a the more important limiting factor to the 
radiative LDOS. 

Owing to the great recent progress in nano-fabrication 
techniques, the scale on which metamaterials can be pat- 
terned is entering the nanometer regime, where nonlocal 
response of the metal becomes important For exam- 
ple, nonlocal response can significantly blueshift the lo- 
calized surface plasmon polariton (SPP) resonance peak 
and modify the field enhancement of a nanoscale plas- 
monic structure^£ ~ 22 i 26 

In this paper, we discuss the effects of nonlocal re- 
sponse on the optical properties of hyperbolic metama- 
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FIG. 1: (Color online) Sketch of a multilayer hyperbolic meta- 
material consisting of periodic dielectric-metal bilayers. The 
dielectric and metal layer thicknesses are ad and a m , respec- 
tively, and their sum equals the period a of the unit cell. 
Corresponding permittivities are td and e m . The red arrow is 
a dipole emitter located in the middle of a dielectric layer. 



terials. It is shown that the nonlocal response gives 
rise a large-wavevector cutoff in the dispersion, inversely 
proportional to the Fermi velocity of the electron gas. 
In fact, the dispersion in hyperbolic media becomes no 
longer strictly hyperbolic. Accordingly, we identify a new 
and fundamental limit on the enhancement of the radia- 
tive emission rates of HMMs. In particular, we show that 
the radiative LDOS does not grow arbitrarily large even 
in the ideal limiting case that all three aforementioned 
parameters 7, a, and D vanish, since the intrinsic non- 
local response turns the 'super-singularity' into a finite 
broadband LDOS enhancement. On a more general level, 
our results illustrate the need, as for metallic nanopar- 
ticle o 25 ' 26 , to take nonlocal response into account in ho- 
mogenization theories, where the goal is to predict the 
effective properties of metamaterials with ever decreas- 
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ing unit cell sizes. 

The paper is organized as follows: In Sec. II, we intro- 
duce the linearized hydrodynamic Drude model within 
the Thomas-Fermi approximation. In Sec. Ill, the un- 
usual dispersion curves of the HMMs and their effective 
parameters are discussed. In Sec. IV, to understand bet- 
ter the HMM dispersion found in Sec. Ill, we discuss the 
SPP supported by a single metal layer. In Sec. V, we in- 
vestigate the LDOS of the HMMs, before discussing our 
results and concluding in Sec. VI. Finally, details of the 
calculations can be found in Appendices A-C. 



II. HYDRODYNAMIC DRUDE MODEL 

We consider a similar multilayer HMM geometry as 
in the recent experiments by Tumkur et al*£, see Fig. [T] 
The unit cell is a sub- wavelength dielectric-metal bilayer, 
which is relatively simple and cheap to fabricate^ and 
allows analytical analysis. For an effective-medium de- 
scription of such a metamaterial, a local-response ap- 
proximation (LRA) is usually employed, i.e. spatial dis- 
persion is neglected. This gives the effective dispersion 
relation 
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hydrodynamic Drude model (HDM) within the Thomas- 
Fermi approximation^ - — In the HDM, the metal sup- 
ports both the usual divergence-free ('transverse') and 
rotation-free ('longitudinal') waves. Above the plasma 
frequency both types of waves can propagate. The 
dispersion fcr(w) of the transverse waves is given by 



follows from £^(fe, oj) 
tions 



while k^iuj) of the longitudinal waves 



0, in terms of the dielectric func- 
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+ iujj — (3 2 k 2 



(2a) 



(2b) 



Here, 7 is the Drude damping, ui p is the plasma fre- 
quency, and the nonlocal parameter (3 is equal to -\/3/5uf 
with vf representing the Fermi-velocity. While e]^ is the 
familiar Drude dielectric function, ej^ depends on up and 
describes nonlocal response. 



III. DISPERSION AND EFFECTIVE 
MATERIAL PARAMETERS 
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Below the plasma frequency cu p , 



where el,° c = a(ad/ed - 
and fen = (fc 2 + fc2)V2 
where e m < in the Drude model of a pure plasma, the 
dielectric tensor elements e zz and e m can have opposite 
signs by a proper choice of the filling factor a m /a. Then 
the dispersion becomes hyperbolic, meaning that an iso- 
frequency contour becomes a hyperbola rather than the 
usual ellipse in the (fe 2 , fey )-plane. The length of this con- 
tour diverges, and so does the radiative LDOS. Although 
we discuss metal-dielectric bilayer structures, we want to 
point out that our theory may also be applied to struc- 
tures where the metal is replaced by other materials with 
a Drude response^ 

In the present paper, we go beyond the LRA, and dis- 
cuss the optical properties of the HMMs in the linearized 



To calculate the exact dispersion equation for the in- 
finitely extended HMM, we employ a transfer-matrix 
method for both transverse and longitudinal waves com- 
bined. Our method is quite similar to the one devel- 
oped by Mochan et al£&, but we corrected the additional 
boundary condition (ABC) that in Ref. HH was employed 
for simplicity. An ABC is required to complement the 
usual Maxwell boundary conditions, and all boundary 
conditions together make the solution to the coupled 
Maxwell and hydrodynamic equations unique. Details 
how to derive the correct ABC and a consistency check 
can be found in Appendix A. 

For arbitrary unit cell size a and metal and dielectric 
filling fractions, we find the exact dispersion relation for 
the infinite HMM to be 
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with 



% = k z a, 9 d = k dz a dl 9 m = k mz a m 61 = fe mz a m , (4) 
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where k 2 dz + k 2 = u 2 e d /c 2 , (fc^J 2 + k 2 = w 2 e^/c 2 , and 
(fcLjz + k 2 _ k 2 with k 2 _ ( w 2 + 11UJ _ w aj )8 a_ This 

dispersion equation looks very similar to the one found in 
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Ref. l28|, but the essential difference is that the parameter 
Wd here is Wd/frf in Ref. [28|. 

As an example we consider a HMM with free-space- Au 
bilayer unit cell with ad = a m = a/2, and we include the 
Au Drude loss. Fig. [5] depicts HMM dispersion curves at 
uj = 0.2w p . Fig. [Ha) shows hyperbolic dispersion in the 
small wavevector regime, whereas Fig. [2jb) zooms out 
and shows strong deviations from hyperbolic dispersion 
for large wavevectors. 

First, in Fig. [21a) we observe three non-coinciding hy- 
perbolic dispersion curves in the small-A: region, one for 
local and two for nonlocal response. This tells us that the 
HDM is not just a local theory with a large-wave vector 
cutoff added, since then the curves for local and for nonlo- 
cal response would have coincided for small wavevectors. 
Furthermore, the two nonlocal hyperbolic curves do not 
coincide, the one for a strongly subwavelength unit cell 
a = Ap and the other for a — > 0. This illustrates that 
the size of the unit cell affects effective-medium proper- 
ties, even in the deep subwavelength limit, which goes 
against common wisdom obtained in the LRA. The rea- 
son for this is that in the HDM the longitudinal wave in 
the metal layer has a large vector k-^(oj) 2tt/X, so that 
typically the condition |&;L(w)|a <C 1 is not satisfied even 
in the deep subwavelength limit. Thus, the longitudinal 
wave can probe the finite size of the unit cell even though 
a <C A, and this gives rise to the periodicity-dependent 
dispersion curve of Fig. HJa). 

Zooming out, Fig. HJb) shows that nonlocal response 
gives rise to closed non- hyperbolic dispersion curves, for 
both considered values of a, in stark contrast to the fa- 
miliar hyperbolic curve in the LRA which is also shown. 
(We still call these media hyperbolic because of their hy- 
perbolic small-wavevector dispersion.) Both and k z 
are bounded on the curve for a — Xp. For smaller values 
of a, we do not expect the hydrodynamic Drude model 
to appl y 17 ' 22 , but as we shall see below it is useful to also 
consider the limit a — > 0. The curve for a — > shows 
a turning point at fcy = k? r In the lossless limit, no ra- 
diative modes exist above fejr, as we explain shortly. The 
wavevector kf, is found to be 



kf, = % oc — . 



(6) 



To analyze the dispersion curves of Fig. [5J we derive 
the effective material parameters of the HMM by a mean- 
field theory that can be applied to many geometries. In 
the limit of vanishing unit-cell size, we obtain the effective 
material parameters 
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where e zz and represent the effective parameters of 
the metamaterials when the metal layer is replaced by 
a free-space layer, with e zz — a(ad/ed + dm) -1 , and 
ae^ = aded + flm- Both nonlocal effective material param- 
eters of Eq. (J7J differ from the corresponding parameters 
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FIG. 2: (Color online) Dispersion curves of the HMM for 
uj = 0.2o;p, on (a) small and (b) large wavevector intervals. 
Red curves for a = Af, green curves for a — > 0, black curves for 
a —¥ in the LRA. The unit cell of the HMM is a free-space- 
Au bilayer with ad = a m = a/2. Material parameters for Au: 
hiop = 8.812eV, h-y = 0.0752eV, and v F = 1.39 x 10 6 m/s. 



for local response. The derivations leading to Eq. are 
presented in Appendix B. Neglecting loss at first, we find 
from Eq. that e" loc has a resonance at fc y = k^ where 
both e" loc and k z diverge. The value of k^ is independent 
of ed (unlike what one would find when using the in- 
correct ABC of Refsj 25 ' 28 ). Increasing k^ beyond fcj?, the 
e" loc changes sign from negative to positive. Since e^ oc is 
always positive, it follows that no mode exists above fcjj\ 
Thus, nonlocal response gives rise to a large- wavenumber 
cutoff at fc|| = k^. With loss, the resonance is smoothed 
out and modes exist also above fejr. However, for fcy — > oo, 
the corresponding k z approaches ioo, which shows that 
such large- wavevector modes are purely evanescent. This 
explains why the dispersion curves in Fig. [5] are closed. 

We stated in Eq. that unlike in the LRA, in the 
HDM the effective parameter ef° c simply equals the 
(positive) permittivity e A zz . This outcome is fixed for 
a — > by the continuity of the normal components of 
the displacement field and the ABC of Eq. (|A1|) with 
£ othcr _ 2 1 7 In particular, the different boundary condi- 
tions explain why the local and nonlocal a —¥ curves 
in Fig. (Ha) exhibit different hyperbolic small-wavevector 
dispersion. 

Above the plasma frequency, the HDM and the LRA 
also exhibit qualitatively different dispersion. In the LRA 



no hyperbolic dispersion exists for frequencies above the 
plasma frequency, not even for small wavevectors, since 
then both ed and e m are positive. By contrast, hyperbolic 
dispersion can exist in the HDM for io > uj p , because 
the effective-medium parameter e^ loc given in Eq. can 
assume negative values above u> p . 



IV. SURFACE PLASMON POLARITON 
SUPPORTED BY A SINGLE METAL LAYER 

In Sec. Ill, it was demonstrated the dispersion curves 
in the LRA and HDM differ significantly. To understand 
this better, here we relate these essential differences to 
the different properties of single metal layers in both the- 
ories, knowing that the bulk modes of the HMM result 
from the coupling of SPPs of neighboring metal layers. 
So we investigate the SPPs supported by a single metal 
layer, first analytically in the quasi-static limit. With 
respect to the magnetic field, the SPPs can be classified 
as even and odd modes. In the HDM, the dispersion 
relations of the even and odd modes are found to be 
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k sp ' k p 



where fc sp represents the SPP wavevector, and k\ z — 
(fc 2 p — fc^) 1 / 2 . In the limit a m — > 0, the dispersion equa- 
tion of the even mode has no solution, but the odd mode 
always has one, even above w p . Its dispersion is such that 
fc sp has k^ as an upper bound in the limit a — ¥ 0. So we 
can now understand that it is this nonlocal 'ceiling' for 
the single-layer SPP wavenumber that leads to a cutoff 
of fc|| for the bulk modes of the metamaterial, as we saw 
in Fig. dJ 

In Fig. [3] we analyze numerically the effect of retarda- 
tion on the SPP dispersion of a single Au layer in free 
space, for local and nonlocal response. With retarda- 
tion, near the light cone also even-mode solutions ex- 
ist. Only for nonlocal response do we find modes above 
w p . Again we find that nonlocal response gives rise to 
a forbidden region k sp > k^ for the odd SPP mode, see 
Fig. By contrast, in Fig. Qfa) for the LRA, both 

even and odd modes have finite-frequency solutions with 
k sp approaching infinity, which leads to the characteristic 
hyperbolic curve of the HMM that extends to infinitely 
large wavevectors. 



V. LOCAL DENSITY OF STATES 

The discussed dramatic modification of the metamate- 
rial dispersion due to nonlocal response will also strongly 
affect the broadband super-singularity known to occur in 
the local-response LDOS, as we shall see. In general, the 
LDOS is proportional to the spontaneous-emission rate 



FIG. 3: (Color online) Dispersion curves of the SPP mode 
supported by a single lossless Au layer with a thickness a m 
in free space, in the (a) local-response aproximation, and (b) 
the hydrodynamic Drude model. Dashed and solid curves 
correspond to even and odd modes, respectively, with red 
curves for a m = 0.1A P , green for a m — 0.01A P , and black 
curves for a m — > oo (single-interface SPP). The gray areas 
are forbidden regions for the SPP modes, with light cones on 
the left. 



averaged over all solid angles, and defined as 

LDOS(r ,w) = -^Tr{lm[G(r ,r , W )]}, (9) 

where G is the dyadic Green function of the medium and 
ro the position of the emitter. The Green function G is 
defined by 

-VxVxG(r,r') + fc 2 J dri e(r, n)G(ri, r') = l<5(r-r'), 

( 10 ) 

where I represents the unit dyad, and e represents the 
dielectric function, which is a position-dependent delta 
function for the local dielectric medium, and a tensorial 
nonlocal operator defined by Eq. ([2]) for the metal. For 
the multilayered HMM, G can be decoupled into separate 
contributions from TM and TE modes. TM modes sup- 
port the hyperbolic dispersion curve, and greatly domi- 
nate the LDOS, so we will neglect the TE contribution 
to the LDOS. 
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If we first neglect loss, then only radiative modes con- 
tribute to the LDOS. For an electric dipole with moment 
H, the contribution to the LDOS of a single radiative 
mode is proportional to \fi- ak(ro)| 2 /|Vk^|, where ak is 
the properly normalized mode function^ In the LRA, 
for the limiting case of a — > 0, the single-mode contribu- 
tion to the LDOS scales linearly in k as fey and k z tend to 
infinity. This results in a diverging radiative LDOS, the 
broadband LDOS supersingularity of hyperbolic media. 

Let us now consider the LDOS in the HDM instead. If 
we again take the limit a — > 0, and let fey tend to and 
k z to infinity, then this time the single-mode contribution 
to the LDOS scales as l/fc z 2 , which we derived using the 
effective parameters of Eq. (0. Radiative modes with 
large wavenumbers are therefore negligibly excited. As a 
main result of this paper, we consequently find that in 
the HDM the radiative LDOS converges to a finite value 
as a — > 0, even though the integration area in k-space 
diverges. We find the numerically exact value and its 
analytical approximation 



LDOS(w) 



6tt 2 ^ 3 



(11) 



where 



n / 2 co- 
de 



i 2 e + (4/4 s )[sm*e-<?rrt] 



sin 



c loc 



/e d 



(12) 

with #o equal to arcsin(e)| OC /e^) for ej° c > and vanish- 
ing otherwise. The derivations leading to Eq. (fTT|) arc 
presented in Appendix C. As illustrated below, Eq. (fTTj) 
entails that nonlocal response leads to a large upper 
bound to the radiative LDOS of the HMM, proportional 
to uj 2 /vp. This exceeds the free-space radiative LDOS 
approximately by c 3 /vp, which is of order 10 7 for most 
metals. 

When taking metallic Drude loss into account, then 
the LDOS has contributions both from radiative modes 
and from nonradiative quenching, the latter due to loss. 
For the limiting case of a — »• 0, we already discussed that 
e" loc tends to see Eq. Q. For large wavevectors fey 
also the other component e"l° c tends to e^ z . Thus, to 
the extent that ed is lossless, the evanescent mode with 
large fey does not contribute to the nonradiative LDOS, 
which therefore stays finite. As a result, the total LDOS 
containing both radiative and nonradiative contributions 
in the HDM converges as a — > 0. In the low-loss case, 
where the radiation LDOS is dominant, Eq. (fTTj) is an 
accurate expression of the total LDOS, as we verify by 
numerically exact simulation below. 

We calculate the LDOS numerically exactly by merg- 
ing two methods: the local-response transfer matrix 
method by Tomas to calculate the Green function of arbi- 
trary multilayer medial, and the aforementioned HDM 
extension of the transfer matrix method^ The details 
can be found in Appendix D. 

Figure^a) depicts the LDOS enhancement, defined as 
the ratio between LDOS in the HMM and in free space, 



as a function of the periodicity a. Clearly, in the HDM 
the LDOS converges to a finite value as a — > 0. This 
proves that both the radiative and nonradiative LDOS 
in the HDM are finite. By contrast, in the LRA the 
LDOS diverges as 1/a 3 , where the nonradiative LDOS 
has a dominant contributioni&2ir— The HDM ceases to 
be valid for a < Af , but the LDOS enhancement value in 
the limit a — > is a useful upper bound. 

We also calculated the LDOS for the lossless case [not 
shown Fig.@Ja)], and we find smaller values for the LDOS 
owing to the missing nonradiative contribution, but the 
same trend for a — > 0. This proves that nonlocal response 
rather than loss is responsible for removing the singu- 
larity of the radiative LDOS. In Fig. @|b) we compare 
the limiting LDOS for the numerically exact method in 
the lossy case with the lossless analytical approximation 
of Eq. (jlll) , when artificially varying the Fermi velocity. 
The value from the exact method is only larger than that 
from Eq. (ITTj) by around 6%. We attribute the small dif- 
ference in LDOS to the nonradiative LDOS due to the 
Drude loss. Thus, the loss acts as a small perturbation 
to the radiative LDOS. 
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FIG. 4: (Color online) (a) LDOS versus the periodicity a at 
the center of the free space layer of the hyperbolic metama- 
terial for uj = 0.2w p . (b) The a —> limiting value of the 
LDOS in the hydrodynamic Drude model as function of vf /c. 
Parameters of the hyperbolic metamaterial as in Fig. [2] 



VI. DISCUSSION AND CONCLUSIONS 

For finite-sized unit cells, small loss gives rises to a 
regular perturbation of the radiative LDOS, both in the 
local and in the nonlocal response theories. However, 
the theories start to differ dramatically in the limit of 
infinitely small unit cells. In particular, in the nonlocal 
hydrodynamic Drude model the small variation of the 
radiative LDOS with small loss is quite different from 
the previously found radiative LDOS scaling with loss in 
the local theory as 7~ 3 / 2 for infinitely small unit cells^ 

The small increase of the total LDOS due to loss is also 
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quite different from spontaneous-emission rates of a point 
emitter inside a homogeneous absorbing medium, where 
the loss induces non-radiative quenching that can dra- 
matically decrease the radiative decay efRciency^^ir— 
In this sense, the nonlocal response regularizes the sin- 
gularity not only of the radiative but also of the nonradia- 
tive LDOS of a lossy HMM. One can interpret this finite 
nonradiative LDOS as due to a nonlocal screening of the 
electron scattering lossi^i In a certain high wave-vector 
region the reverse can also occur, namely the enhance- 
ment of the nonradiative LDOS by the nonlocal response, 
when not only taking Drude loss into account, as we do 
here, but also electron-hole pair absorption. By neglect- 
ing any dielectric response of the metal apart from the 
(hydrodynamic) Drude response, we underestimate the 
nonradiative LDOS of real metals. However, the impor- 
tant conclusion that the nonlocal response removes the 
singularity of nonradiative LDOS is still valid<24 

In conclusion, we have shown that the hydrodynamic 
Drude model gives closed non-hyperbolic dispersion re- 
lations for hyperbolic metamaterials, with a fundamen- 
tal wavevector cutoff cx oj/vp. These effective disper- 
sion relations have hyperbolic limits for small wavevec- 
tors, but the precise hyperbola depends on the subwave- 
lcngth size of the unit cell, contrary to consensus based 
on the local-response approximation. We find that the 
hydrodynamic model regularizes the broadband super- 
singularity of the radiative LDOS, and provides a large 
physical upper bound proportional to uj 2 /vp. In practice, 
considering the finite values of a and D, i.e., the finite 
sizes of the unit-cell and the emitter, we usually have 
1/a < 1/D < uj/v-p. This indicates that the size effects 
have a dominant role in limiting the LDOS enhancement. 
Thus, under an upper bound set up by the nonlocal re- 
sponse, hyperbolic metamaterials have a plenty of room 
for improvement in boosting light-matter interactions by 
decreasing the sizes of the unit-cell and the emitter. 
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Appendix A: Boundary Conditions 

Since the hydrodynamic dynamics allows the exci- 
tation of longitudinal waves, the unambiguous solu- 
tion of the nonlocal-response dynamics requires ad- 
ditional boundary conditions (ABCs), complementing 
the Maxwell boundary conditions. As is well known, 
the Maxwell boundary conditions are a consequence 
of Maxwell's equations themselves, in the sense that 
the derivation of the boundary conditions only involves 
Maxwell's equations plus mathematics (the Gauss and 
Stokes theorems). Quite analogously, ABC's are not a 



matter of choice but can be derived from the (linearized) 
hydrodynamic equations, at least for a given equilibrium 
free-electron density profile noi 17 ' 35 When assuming a 
simple zero-to- nonzero step profile of hq at the dielectric- 
metal interfaces, this unambiguously leads to one and 
only one required ABC, namely the continuity of the nor- 
mal component of the free-electron current Ji 17 i 35 

Let us now write the relative permittivity of the di- 
electric medium as £d, and the dielectric response of the 
metal as e m (ui). We assume that e m (ui) is given by the 
sum of a nonlocal hydrodynamic Drude free-electron re- 
sponse plus e^ her (w), the latter describing the remain- 
ing dielectric response of the metal. Since one of the 
Maxwell boundary conditions is the conservation of the 
normal component of the displacement field, the ABC is 
equivalent to the condition 



£ ther E m • n = e d E d • n, 



(Al) 



where E m .d represent the electric fields in the metal and 
dielectric, respectively, and h is the unit vector normal 
to the boundary. From Eq. (|A1[) , we see that the normal 
electric field is discontinuous across the boundary when 
£ other £rf j um p m the electric field occurs due to 
the surface charge produced by polarization of the bound 
electrons both in the dielectric and in the metal. 

We discuss the ABC in some detail, because Mochan 
et a/. 28 , whose pioneering transfer matrix method we em- 
ploy here, and also recently Ciraci et al^- used instead 
the continuity of the normal component of the electric 
field as the ABC, 



E m • n = E d • n, 



(A2) 



or equivalently the continuity of the normal component 
of the displacement current. There is no derivation of the 
latter ABC in Refs. 25,28 . It happens to be only correct, 
in agreement with Eq. (jAll) . if e^ 1101 ' = £d, for example 
in case the dielectric is vacuum (ej = 1) and the metal is 
a pure Drude metal (e^ 1101 ' = 1)- Mochan et al£2- applied 
their ABC for simplicity and write that they thereby ig- 
nore the discontinuity of the electric field, due to the 
accumulation at the surface of bound charges. Our main 
point is here that without additional complication the 
correct ABC can be implemented, and that many phys- 
ical predictions of the hydrodynamic Drude model are 
sensitive to implementing the ABC correctly. 

To understand the ABC physically, recall that in the 
HDM the dynamics of the free electrons is described by 
the equation of motion 



9v _ 



Vp. 



'dog 



e(E + vxB), (A3) 



where pdeg is the pressure from the ground state en- 
ergy of the degenerate quantum Fermi gas, and n is 
the free-electron density. The pressure force —'Vpdeg/n 
oc -Vn/n drives the free electrons diffusing from the 
high density region to the low density region. It is this 
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force that prevents the free-electron charge from accu- 
mulating on the boundary surface, since the existence of 
a free-electron surface charge would cause an infinitely 
large pressure force, which is unphysical. The nonexis- 
tence of the surface free-electron charge indicates that 
the free-electron current should be continuous across the 
boundary, as the ABC (|A1|) indeed describes. By con- 
trast, in the ABC of Eq. (|A2I) . there exists no surface 
charges at all. This indicates that the pressure force 
somehow smears out not only the surface free-electron 
charge in the metal but also the surface polarization 
charges in both the metal and the dielectric. However, 
one cannot expect the smearing out of the surface po- 
larization charges in the HDM, since the pressure force 
only acts on the free electrons in the metal. In this sense, 
the ABC of Eq. (|A2I) is not consistent with the assumed 
dynamics and thus not physically sound. 

There is another perhaps simpler argument, a consis- 
tency check that confirms that the ABC of Eq. (|A2|) is 
more problematic. Assume there is a thin free-space 
layer with subwavelength thickness S between the non- 
local metal and the local dielectric medium. At the 
boundary between the metal and free space, the ABC 
of Eq. (|A2[) gives E m ■ n = Ef • n, where Ef represents 
the electric field in the free-space layer. At the boundary 
between free space and the dielectric medium, we have 
Ef • h = fidEd • n by the standard Maxwell boundary 
condition of the continuity of the normal component of 
the displacement field. In the limit of an infinitely thin 
free-space middle layer (S — > 0), the three- layer system 
essentially becomes the two-layer system where the metal 
and the dielectric medium touch, and for which we find 
E m • n = e^Ed ■ n by combining the previous two iden- 
tities. However, this contradicts with Eq. (| A2|) for the 
metal-dielectric interface. Thus, the ABC of Eq. (IA2|) 
can not be applied consistently. For the ABC of Eq. (|A1[) , 
we obtain instead consistent results when following the 
above thin-layer argument. 



Appendix B: Effective material parameters of 
hyperbolic metamaterial 

When the unit cell has a thickness a that is much 
smaller than an optical wavelength Ao, then the optical 
properties of such an infinite multilayer structure can be 
macroscopically described by a diagonal effective dielec- 
tric tensor e = diagfey , ey , e zz ] with tensor components 



■'■■!:) I 



(E, 



(Bl) 



and where (. . .) denotes spatial averaging over a unit-cell. 

The unit cell can be chosen symmetric, identical for 
left- and right-traveling waves. Consider a unit-cell po- 
sitioned at —a/2 < z < a/2, with the metal layer at 
— a m /2 < z < a m /2, which is symmetric in z — 0. The 
total fields in such a unit cell are generated by waves 
incident both from the left ("1") and from the right 



"r", and the average fields can be split into two terms, 
(E) = (E)i + (E) r . However, by symmetry of the unit cell 
it follows that e = (D)/(E) = <D)i/(E)i = (D) r /(E) r . 
To obtain the effective material parameters, we can sim- 
ply replace the average fields in the periodic structure by 
the average fields in a single unit cell. 

Before spatially averaging the fields, we first need to 
find them as solutions of Maxwell's equations. We focus 
solely on TM-polarized waves since the hyperbolic disper- 
sion occurs for those waves only. Since kaa m <C 1, we can 
make the quasi-static approximation, where E = — V0. 
Consider an incident electric field with the electric po- 
tential 4> = eyLp{ik\\X — k^z). The electric potential in the 
whole system can then be written as 

4>\ = exp(?fc||2;) [exp(— fcyz) + r exp^k^z)] , 

4>2 = exTp(ik l{ x) [Ai exp(— fc|,z) + A 2 exp(fcnz)] , 

02 = exp(ifcna;) [Bx exp(-fc Lz z) + B 2 exp(fc Lz z)] , 

03 = t exp(ifcj|2;) exp(— k^z), 

(B2) 



where khz — — ^l- By matching boundary con- 
ditions at the two metal-dielectric interfaces, the above 
equations can be solved. After obtaining the field distri- 
butions, we can average the fields from —a/2 < z < 
a/2 to obtain the effective material parameters using 

Eq. 4ED- 

First consider the case in which the metal layers are 
much thicker than the wavelength of the longitudinal 
waves (/cL Z a m ^ 1), but where the unit cell is much 
thinner than an optical wavelength. Here we expect an 
effective (homogenized) description to apply and nonlo- 
cal response to be negligible. Following the above de- 
scribed scheme, the effective material parameters, to the 
zeroth-order in the small parameters k^a and l/(fcL Z a m ), 
are found to be 



doc 



1 



id 

1(1 



e,! 00 — /d«d + / m e„ 



(B3) 



in terms of the filling factors /d = a^/a and / m = a m /a. 
Indeed, the effective material parameters are just as 
what one would find in the local response approximation 
(LRA). 

Second, we consider the limiting case of a m — > with 
khzO-m <C 1, where the nonlocal response is extremely 
strong. As before, we keep the filling fractions / m ,d con- 
stant when taking the limit. The effective material pa- 
rameters, to zeroth order in both k«a and fcL Z a m , become 



L.2 Joe/ d _ 1,2 T 
nloc d nloc _ ri ft 'L c ll / c ll n '|| c m 



^zz ^zzi 



h 2 — h 2 f T 



(B4) 



Eq. (|B4|) characterizes how the nonlocal response can 
modify the effective material parameters to the largest 
extent. 



8 



Appendix C: Limiting LDOS of hyperbolic 
metamaterials 



Here we provide the calculation details of the LDOS in 
the hydrodynamic Drude model in the limit of infinitely 
small unit cells (a — >• 0). We employ the effective material 
parameters derived in Eq. (|B4I) . Since the by far domi- 
nant contribution to the LDOS stems from TM waves, we 
will neglect the TE contribution. In fc-space, the diagonal 
components of G, in an effective medium with material 
parameters expressed in Eq. (|B4|) for TM polarization, 



are found to be 



J TM, 



jU2^nloc 1*2 L2 /^iiloc U2 / c nlo 

ft n e n K o / e -zz ^z/ e n 
1 I - kl/{klef° c ) 

^nloc L.2 L.2 /^nloc U2 /^nloc ' 

t zz "0 HI / c zz "'z/ t || 



for j = x, y, 
(CI) 



When inserting these diagonal components for the Green 
tensor into expression Eq. (O for the LDOS, we obtain 



lim LDOS 

a->0 



2fcn 1 



3ttc (2 

6n 2 c 
fc 



l^Im J d 3 k G 



k 

TM, jj 



Re 



67T 2 C 



6tt 2 /3 : 



Re 



dfci. 



oo L.2 

dfc 



^11 fe* 



fc 3 i_ ( ;.oc /e J= 



L.2^nloe jU2 £ nloc /^iiloc 

fc || fc || / t zz 

l-ef°°/e?r 



kn I ^nloc /^iiloc 



(C2) 



where ?y is expressed in Eq. (|12[) . In the derivation, 
we neglected losses in the metal. To arrive at the sec- 
ond line of Eq. (|C2I) , we use the principal- value identity 

lim -At = P- ± iirS(x). The final identity then follows 

A4.0 x±ta x 

immediately by inserting the a — > limiting expressions 
for ef° c and e" loc given in Eq. (|B4t . 



Appendix D: Green Function of hyperbolic 
metamaterial 



sents the Green's functions owing to the scattering be- 
tween central layer and the left and right semi-infinite 
HMM. In the plane wave basis, Gd is expressed as^S, 



G d (r, r ) 



5{z - z ) 
fed 2 



d 2 k y exp [ik„ • (r,| -r „)] 



i f , 2 i [ e TE e TE + eTMeTM 1 ] 

exp[ik|, • (r,| - r ||) + ik z \z - z \], (D2) 



Consider an emitter positioned in the dielectric layer of 
the HMM. The HMM can be divided into three regions: 
(1) the central dielectric layer where the emitter is lo- 
cated; (2) the left semi-infinite HMM; (3) the right semi- 
infinite HMM. The distance between the emitter and the 
left (right) boundary of the dielectric layer is zi (z r ). The 
Green function G in the central layer could be separated 
into two terms 



G(r, r ) = G d (r, r ) + G s (r, r ), 



(Dl) 



where Gd represents the Green's function for the emitter 
in the homogenous dielectric medium, while G s repre- 



with 



e TE = jr x z 



h 

fey 

ku i h z z 



"TM 



kd 



x e T E, 



(D3) 



where k y = k x x + k y y, kd = oj^/ed/c, fcj + jfe? = fc 2 , and 
e^ M correspond to z > zq and z < zq, respectively. The 
terms containing exE and exM represent TE and TM 
waves, respectively. 

The scattering part G s of the Green function is ex- 
pressed as 
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i f 1 

G s (r,r ) = ^2 / dk|,— exp[ik,| • (r„ -r |,)] 

X { [^*TE e TE e TE + r ™ e TM e TM + r TE e TEe T E + r ™ e ™ e TIVl] ex P i ik z( z ~ z o)} 

+ [ r TE e TEe T E + r ™ e ™ e TM ~^TE e TE e TE + ^™ e ™ e TM] ex P HM Z ~ Z )]} , (D4) 



where r ±=t is the reflection coefficient, in which the left 

TE,TM 1 

superscript "±" represents the scattering wave in the ±z 
direction, and the right superscript "±" represents the 
incident wave in the ±z direction. The r ±:t are found 

TE.TM 

to be 



^TE TM G x p(2ifc z O.(j) 



' TE,TM 



' TE,TM 



' TE,TM 



' TE,TM 



1 _ ^te,tm exp(2ifc z a d ) ' 

-RTE,TMexp(2ifc z 2: r ) 



1 _ ^te,tm exp(2ifc z a d ) 

^TE,TMexp(2zfc z Z;) 

1 - i?2 e tm ex p(2ifc z a d ) 



(D5) 



where -Rte,tm denote the reflection between the dielec- 
tric medium and the semi-infinite HMM for TE and TM 
waves, respectively, z r (zi) represents the distance be- 
tween the point emitter and the right (left) semi-infinite 
HMM. -Rte.tm can be calculated by the transfer matrix 
method as demonstrated in Ref. [23. 
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